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In the quest for quantum spin liquids in three spatial dimensions (3D), we study the half-filled 
Hubbard model on the simple cubic lattice with hopping processes up to third neighbors. Employing 
the variational cluster approach (VGA), we determine the zero-temperature phase diagram: In 
addition to a paramagnetic metal at small interaction strength U and various antiferromagnetic 
insulators at large U, we find an intermediate-t/ antiferromagnetic metal. Most interestingly, we 
also identify a non-magnetic insulating region, extending from intermediate to strong U. Using 
VGA results in the large-U limit, we establish the phase diagram of the corresponding J1-J2-J3 
Heisenberg model. This is qualitatively confirmed - including the non-magnetic region - using spin- 
wave theory. Further analysis reveals a striking similarity to the behavior of the J 1 -J 2 square-lattice 
Heisenberg model, suggesting that the non-magnetic region may host a 3D spin-liquid phase. 


Introduction. — Quantum spin liquids (QSLs), charac¬ 
terized by the absence of long-range order down to the 
lowest temperatures, are amongst the most fascinating 
states of matter mm- Their elementary excitations of¬ 
ten display fractionalization, resulting in unconventional 
quasiparticles with exotic properties. QSL states typi¬ 
cally emerge from a combination of frustration, i.e., the 
competition of different exchange couplings, and quan¬ 
tum fluctuations, both acting to suppress symmetry¬ 
breaking order. Given that fluctuation effects in lat¬ 
tice systems weaken with increasing coordination num¬ 
ber, the majority of works in the field have focussed on 
systems in two dimensions (2D), while 3D spin liquids are 
scarce because they are harder to stabilize. A notable ex¬ 
ception are spin-ice pyrochlores which, however, are clas¬ 
sical spin liquids [5]. Concerning 3D QSLs, some of the 
few promising experimental candidates are the hyperka- 
gome compound Na4lr308 [ 1 ] and the pyrochlore system 
Pr2lr2O7 0, but also Tb2Ti207 i, Yb2Ti207l7], and 
FeSc2S4 [8] might fall into this category. 

Theoretical work is mainly limited to a few exactly sol¬ 
uble spin models [SHU], quantum dimer models mmi, 
and spin-liquid states predicted using parton construc¬ 
tions [HHIll ■ In contrast, results from numerical sim¬ 
ulations, proven to be extremely powerful in the 2D 
case, are largely missing for 3D. The main reason 
is computational complexity, which severely limits the 
application of exact-diagonalization and density-matrix 
renormalization-group methods. While quantum Monte 
Carlo techniques are available in principle [19] , models of 
frustrated electrons typically suffer from the sign prob¬ 
lem. Relevant models in 3D have therefore been stud¬ 
ied using approximate methods such as spin-rotation- 
invariant Green’s functions |20j , dynamical mean-field 
theory (DMFT) pT][2^ and its cluster extensions |23l I24j 
as well as cluster perturbation theory [25]. However, ob¬ 
taining reliable results concerning magnetic ground states 
has proven notoriously difficult. 


In this Rapid Communication, we partially fill this gap 
by generalizing the VGA [55] to 3D and applying it to the 
Hubbard model with longer-ranged hoppings on the cubic 
lattice. We are able to compute zero-temperature phase 
diagrams and single-particle spectra for arbitrary degree 
of frustration and for arbitrary interaction strength. We 
determine the location of the metal-to-insulator transi¬ 
tion (MIT) as well as the boundaries of various types of 
collinear antiferromagnetic order. Most remarkably, we 
identify an insulating region extending from intermedi¬ 
ate to strong interactions which is devoid of magnetic 
order, rendering it a spin-liquid candidate. We inves¬ 
tigate the large-G limit using VGA and compare this to 
spin-wave theory for the corresponding frustrated Heisen¬ 
berg model, with good qualitative agreement. Based on 
striking parallels to the behavior of the frustrated square- 
lattice Heisenberg model we suggest that the cubic-lattice 
non-magnetic insulating region may host both a 3D QSL 
and a valence-bond solid. 

Model. —We study the Hubbard model with hoppings 
up to third neighbors on the simple cubic lattice, with 

T~L — ~ ~ ^2 ~ ^3 

iu) <(U» {{(u))) 

i 

in standard notation. The chemical potential is chosen 
such that half-filling is ensured. Throughout the paper, 
we assume ti = 1 and 0 < ^2,3 < I- For ^2 = 0 the 
spectrum is particle-hole symmetric. Non-zero t 2 intro¬ 
duces frustration in the sense of destructive interference 
of hopping processes; this is also apparent in the large-17 
limit where t 2 generates a frustrating exchange coupling. 

Variational Cluster Approximation. —VGA is a quan¬ 
tum cluster approach, designed to compute single¬ 
particle properties of interacting many-body systems [26] . 
One first solves a small cluster exactly and derives the 
corresponding full Green’s function. Using the frame- 


2 


(a) is = 0 


AF(7t,7t,n) 

AF(0,7r,7t) 

insulator 

insulator 


AF(o,ii,7i) metal 



paramagnetic metal 


0 0.2 0.4 0.6 0.8 1 

hlh 


(b) ts/ti = 0.5 




NN 

- 

AF(7t,7r,n) 


■ AF{o,7t,n) 

insulator 


insulator 



AF(0,7t,7r) 

metal 

.. paramagnetic 

metal 


0 0.2 0.4 0.6 0.8 1 

h/h 



FIG. 1: Phase diagrams of the cubic-lattice Hubbard model, plotted as function of UjW and 12 /1\ at fixed ts = 0 (a) and 
ts/ti = 0.5 (b), with W being the bandwidth of the dispersion at t/ = 0. The small-t/ paramagnetic metal is destroyed in 
favor of antiferromagnetic insulating phases, for labels see text. In (b) we also find a Mott-insulating regime without magnetic 
order (white region) which persists up to large U. Finally, the shaded region in (a) and (b) is an antiferromagnetic metal, (c) 
Illustration of the three magnetic orders. 


work of self-energy functional theory |27j , one obtains the 
full Green’s function of an infinite lattice which is covered 
by the clusters, and the individual clusters are coupled 
by hopping terms only, but in a self-consistent, varia¬ 
tional scheme. The latter step represents a significant 
approximation to the full many-body problem, while the 
method still includes spatial quantum correlations. Em¬ 
bedded into a grand-canonical ensemble, the stationary 
configuration with lowest potential VL is found by varying 
the chemical potential as well as all single-particle param¬ 
eters of the reference cluster. Magnetic instabilities can 
be investigated by means of Weiss fields. VGA was first 
introduced for ID systems [55] and later generalized to 
2D |5^. In recent work, we have demonstrated that VGA 
is capable of detecting |29| or rejecting |30| non-magnetic 
insulator phases which might correspond to spin liquids. 
For further VGA details we refer to Refs.1551 and 1501 

Here we generalize and apply VGA to 3D for the first 
time. We use clusters consisting of 2^ = 8 lattice sites 
(cube) or 3 X 2^ = 12 lattice sites (rectangular cuboid) to 
compute the full Green’s function of the Hubbard model 
0 - In fact, most data were obtained using the cubic 8 - 
site cluster, and a few selected data points were double- 
checked with the 12 -site cluster: the 12 -site-cluster re¬ 
sults are shown as orange diamonds in Fig.[^(a) and in 
Fig.i deviations turned out to be negligible. As shown 
earlier |30|, it is crucial to vary all the single-particle pa¬ 
rameters, i.e., not only the chemical potential but also 
H, 2,37 and we follow this protocol here. 

Results .—Using VGA we compute the Green’s function 
for each parameter triple [t^ltx^ ta/ti, U/ti) from which 
we extract the single-particle spectral function A(fc,a;). 
Its momentum integral readily allows to distinguish be¬ 
tween metallic and insulating phases. In addition, we 
used Weiss fields to test the presence of the antiferro¬ 
magnetic (AF) orders known from the classical J 1 -J 2 - 


J 3 Heisenberg model: AF(o,^,^), and AF(o,o,, 7 ) 

where the subscript denotes the ordering vector Q [see 
Fig.[^(c)]. These orders are also referred to as G-type, 
C-type, and A-type AF, respectively |3T| . Note that the 
latter two possess a discrete broken Z 3 symmetry in ad¬ 
dition to the broken SU(2) symmetry. Notably, all these 
magnetic orders are collinear and display a two-site mag¬ 
netic unit cell. 

Two representative zero-temperature phase diagrams 
are shown in Figs.[^(a,b). At small U we find the ex¬ 
pected paramagnetic metal, descending from the non¬ 
interacting limit, with the exception of the case t 2 = 0 
where the AF^.^. ,r, 7 r) insulator persists down to infinitesi¬ 
mal values of U due to perfect nesting. While t 2 ^ 0 de¬ 
stroys the perfect nesting and thus stabilizes the small-U 
metal, finite counteracts the effect of t 2 which is re¬ 
flected in an effective reduction of the metallic phase [see 
Fig.[^(b) for t^/ti = 0.5]. At large U we find AF insu¬ 
lating phases for most parameter values, with AF( 7 r, 7 r,.n-) 
(AF(o_.n.^ 7 r)) being realized at small (large) ^ 2 /^ 1 , respec¬ 
tively. If both t 2 and ts are large there is also an AF(q 
insulator. While the MIT coincides with the onset of 
magnetism at small and moderate <2 (and is thus pre¬ 
sumably of first order), for large ^2 we observe instead 
that the MIT is preceded by a magnetic transition: The 
resulting intermediate AF(o_,n._^) metal displays a Fermi 
surface with electron and hole pockets which shrink upon 
further increasing U and disappear at the MIT. 

The most remarkable finding is a narrow insulating re¬ 
gion between AF(,n. ,r, 7 r) and AF(o, 7 r,- 7 r) phases where mag¬ 
netic order is absent. We found this non-magnetic in¬ 
sulating (NMI) region for values of the third-neighbor 
hopping 0 < t^/ti < 0.6. It emerges near the point 
where the paramagnetic metal and two magnetic insu¬ 
lators meet and persists up to the strong-coupling limit 
where it can be associated with a quantum-disordered 
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FIG. 2: Intensity plot of the single-particle spectral function A{k,uj) along the momentum-space path T-X-M-R-T with 
r = (0, 0, 0), X = (tt, 0, 0), M — (tt, tt, 0), and R = (tt, tt, tt). (a) AF(^ o) at parameters = t 2 = 1, ta = 0, and 1/ = IF = 24. 
(b) AF(.^_o,o) for = fa = ts = 1 and C/ = IF = 36. (c) NMI for fa/li = 0.73, ta/ti = 0.5, and U = W = 23.68. For 
comparison, the non-interacting dispersion is shown (solid, red); in (a) and (b), also the Q-shifted (backfolded) non-interacting 
dispersion is shown (dashed, red). 


regime of local moments. Its nature will be discussed 
below. 

Sample results for the single-particle spectrum A(fc, oj), 
measurable in angle-resolved photoemission spectroscopy 
(ARPES), are shown in Fig.l^ AF(^_.,r,o) in panel a), 
AF(,r,o,o) in panel b), and NMI in panel c). The interac¬ 
tion strength U is chosen such that U/W = 1 in all cases. 
In addition to A{k,ui) we also show the non-interacting 
dispersion for comparison. While the data in Figs. [|;a,b) 
correspond to states deep in the AF phase, the resulting 
spectrum is more complicated than a simple two-band 
signal that would be obtained from the unit-cell doubling 
of the band dispersion, indicating Hubbard-type interac¬ 
tion effects. We note that our A(fc, w) data do not display 
true lifetime broadening since we perform VGA calcula¬ 
tion without bath sites. Nonetheless, one can observe 
some precursor of broad bands, i.e., dense sequences of 
levels, which we identify as lower and upper Hubbard 
bands. A remarkable feature in Fig.j^c), corresponding 
to the NMI region, is the weak momentum dependence 
indicating the high degree of frustration - this is already 
visible in the non-interacting dispersion. 

Strong-coupling limit .—In the limit of large on-site re¬ 
pulsion U we expect charge fluctuations to become ir¬ 
relevant. The remaining spin degrees of freedom of the 
Mott insulator are conveniently described using an effec¬ 
tive Heisenberg Hamiltonian which can be obtained per- 
turbatively in t/U. The present Hubbard model features 
real spin-independent hopping, such that every hopping 
process yields an AF Heisenberg coupling Ji = Atf/U. 
Thus the Hamiltonian Q simply reduces to the J 1 -J 2 - 
J 3 spin-1/2 Heisenberg model. Consequently, by fixing U 
to large enough values (here U/ti = 100), our Hubbard- 
model calculation within VGA will yield the J 1 -J 2 -J 3 
quantum phase diagram which has been derived here for 
the first time [see Fig.[^. Since the corresponding phase 
diagram at U/ti = 30 shows only minimal deviations 
compared to Fig.[^ we conclude that we are deep in the 


Mott phase and the Heisenberg description is justified. 
In Fig. 1^ we have exemplarily shown the grand potential 
as a function of J 3 /J 1 at fixed J 2 /J 1 = 0.55 for different 
magnetic solutions (AF(.n._,n. ,r), AF(o_ 7 r, 7 r)i and AF(o,o, 7 r)) 
and the paramagnetic solution (NMI). The stable solu¬ 
tion with lowest value H corresponds to the ground state 
as shown in the phase diagram Fig.|^ 

While there is no alternative method available to test 
our 3D VGA results at intermediate U in the frustrated 
regime, the large-!/ limit lends itself to comparisons with 
approximate solutions of the Heisenberg model. The sta¬ 
bility of the AF phases is most naturally tested using 
linear spin-wave theory (LSWT). This has been applied 
to the J 1 -J 2 model on the cubic-lattice in Ref.Ell and we 
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FIG. 3: Phase diagram of the J 1 -J 2 -J 3 cubic-lattice Heisen¬ 
berg model, obtained by VGA calcnlations for the Hubbard 
model at U/t\ = 100 where J 2 , 3 /Ji = for labels see 

text. Inset: Phase diagram as obtained within linear spin- 
wave theory. The purple dot within the NMI region indicates 
a stable solution of a columnar VBS, see text. 
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FIG. 4: Grand potential f2 as a function of J 3 /J 1 at fixed 
J 2 /J\ = 0.55 for U/t-i = 100. Shown are the energies for 
AF(^ ,r ,r), AF(o,^,.n-), AF(o,o, 7 i-) phase and the non-magnetic 
insulator regime. Ending of lines indicates the disappearance 
of stationary points. 

extend this to include J 3 here. LSWT qualitatively con¬ 
firms the VGA phase diagram Fig.[^ where the LSWT 
result is shown as inset. Most importantly, it confirms 
the existence of the non-magnetic region as found within 
VGA. While the predictive power of LSWT may be ques¬ 
tioned in principle [33] , we point out remarkable similar¬ 
ities to the well-studied case of the J 1 -J 2 niodel on the 
square lattice [M] . Here, LSWT correctly predicts [3S] 
the non-magnetic region around JijJ\ = 0.5 which has 
later been confirmed numerically (see e.g. Refs. [3611401 and 
references therein). Moreover, a quantitative comparison 
of the phase boundaries obtained by LSWT and the 
recent numerical analysis [5^1140] shows that LSWT pre¬ 
dicts a non-magnetic region which is too narrow and also 
shifted towards the AF phase with Q = ( 7 r, 7 r); remark¬ 
ably, that we find the same trend on the cubic lattice 
when comparing LSWT and VGA data. Together, this 
indicates that LSWT is appropriate in benchmarking the 
large-G VGA results [55] . 

The shape of the phase boundaries in Fig. can be ra¬ 
tionalized as follows: Upon increasing J 3 /J 1 the balance 
between the AF(.n. ,^. 71 -) and AF(o_.n._ 7 r) phases is shifted to 
larger J 2 /J 1 - this is because J 3 > 0 does not compete 
with Ji > 0 and thus stabilizes the AF(,r, 7 r, 7 i-) phase rel¬ 
ative to AF(Q ,r,,r). The presence of NMI is related to 
quantum fluctuations: While the competition of J\ and 
J 2 is apparently insufficient to destroy long-range order 
(unlike in the square lattice), the additional fluctuation 
processes arising from J 3 succeed in suppressing order. 

We finally discuss the nature of the NMI region 
which cannot be accessed by the present VGA calcu¬ 
lations. Prominent candidate states are valence-bond 
solids (VBS) with broken translational symmetry and 
various types of spin liquids m- In contrast, long-wave¬ 
length non-collinear spirals (not captured by VGA) ap¬ 
pear unlikely because spirals are absent in the classical 
phase diagram, and quantum fluctuations typically prefer 
collinear instead of non-collinear states. Here we investi¬ 
gate the stability of a columnar VBS, with a two-site unit 


cell and ordering vector (7r,0,0), using bond-operator 
theory |42| . In the harmonic approximation [43] we find 
that the columnar VBS is marginally stable, i.e., it ex¬ 
ists only at {J 2 /Ji, J 3 /Ji) = (0.5,0.25) where it is gap¬ 
less - this is the point where the three magnetic phases 
meet in the classical limit. Again, this reveals a remark¬ 
able similarity to the J 1 -J 2 square-lattice model 
where - within harmonic approximation - a columnar 
VBS state is only stable at the classical transition point 
J2/J1 = 0.5 [15]. It is further known that inclusion of 
interactions stabilizes a gapped VBS state in a finite win¬ 
dow of parameters O 05] , and we expect this to apply 
to the cubic-lattice case as well. While numerical studies 
in the past suggested either a VBS phase or a spin liquid 
to be realized |57H5D] , very recent numerical studies [55] 
of the J1-J2 square-lattice model indicate that its non¬ 
magnetic region is split into a VBS phase (at larger J 2 
in vicinity to the collinear state) and a spin-liquid phase 
(at smaller J 2 in vicinity to the Neel state). Based on 
the striking similarities between the square-lattice results 
and our findings for the cubic lattice we speculate that at 
least a part of the NMI region in Fig. [ghosts a 3D QSL. 
Note that these similarities between 2D and 3D are not 
accidental but rather generic: at the classical transition 
points (J 2 /J 1 , Ja/^i) = (0.5,0.25) (3D) and J2/J1 = 0.5 
(2D) both Hamiltonians can be expressed as the sum over 
operators of elementary cubes (3D) or squares (2D), 
respectively. 

Eventually we comment on the relevance of the pu¬ 
tative spin liquid phase for real materials. To stabi¬ 
lize the non-magnetic insulating phase rather large ra¬ 
tios of J2/J1 and J3/J1 are required. This is, how¬ 
ever, not an obstacle. The search for the spin-liquid 
phase in the analog square lattice case stimulated efforts 
in crystal growth. As a result, the perovskite PbVOs 
features J2/J1 ~ 0.2 ... 0.4 [471 |48|. The compound 
BaGdV 0 (P 04)2 even reaches J2/J1 ~ —0.9 (here Ji is 
ferromagnetic) 05]; and in Li 2 V 0 (Si,Ge )04 couplings as 
large as J2/J1 ^ 12 have been found [501151] . In the light 
of these facts the considered magnetic couplings for the 
cubic lattice might be identified in future experiments. 

Conclusion and Outlook. —We have applied VGA to 
the half-filled frustrated Hubbard model on the simple 
cubic lattice and computed zero-temperature phase dia¬ 
grams. In addition to the weak-coupling paramagnetic 
metals and antiferromagnetic insulators at strong cou¬ 
pling, we have found an antiferromagnetic metal at in¬ 
termediate coupling and - most importantly - an ex¬ 
tended non-magnetic, i.e., quantum-disordered, insulat¬ 
ing region. Using the VGA data at U/ti = 100 we have 
established the phase diagram of the J1-J2-J3 spin-1/2 
Heisenberg model, whose structure could be confirmed 
using linear spin-wave theory. The unexpected but strik¬ 
ing similarity to the J1-J2 square-lattice model suggests 
that the non-magnetic region hosts both a VBS [55] and 
a 3D QSL. 
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Our results raise several issues which will be subject 
of future work: (i) Opposite-sign hopping processes, e.g., 
^ 2/^1 < 0, as realized in cuprate superconductors, may 
modify the physics at intermediate U. (ii) Without per¬ 
fect nesting, superconductivity may occur at half hlling 
at intermediate U. For the cubic lattice this can be tested 
using VGA. (iii) The physics of the Hubbard model Q 
away from half filling is unexplored, but may be tackled 
using VGA as well, (iv) Given that this work establishes 
VGA as a suitable method for 3D Hubbard models, its 
application to other 3D lattices such as hexagonal, py- 
rochlore, and hyperkagome promises to unravel further 
novel physics. 
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